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We introduce a new approach for the numerical pricing of American options. The main idea 
is to choose a finite number of suitable excessive functions (randomly) and to find the smallest 
majorant of the gain function in the span of these functions. The resulting problem is a linear 
semi-infinite programming problem, that can be solved using standard algorithms. This leads 
to good upper bounds for the original problem. For our algorithms no discretization of space 
^ and time and no simulation is necessary. Furthermore it is applicable even for high-dimensional 

I problems. The algorithm provides an approximation of the value not only for one starting point, 

but for the complete value function on the continuation set, so that the optimal exercise region 
CS| and e.g. the Greeks can be calculated. We apply the algorithm to (one- and) multidimensional 

diffusions and to Levy processes, and show it to be fast and accurate. 

1—5 

^\ Key Words: American options, optimal stopping, excessive functions, upper bounds, semi- 

infinite linear programming 

u 

^ 1 Introduction 

I 

, Pricing American type options on multiple assets is a challenging task in mathematical finance and 

is important both for theory and applications. The problem to be solved is an optimal stopping 
fSJ problem. These problems play an important role in many other fields of applied probability, too. 

^ They also appear, for example, in mathematical statistics and portfolio optimization. Although a 

^2 general theory is well developed (cf. e.g. [PS06| ). the value and the optimal strategy in optimal 

stopping problems cannot be found explicitly in most situations of interest. Many approaches 
have been proposed for a numerical solution of optimal stopping problems in the last years. 
For pricing the standard American options in the Black-Scholes market with one underlying the 
most prominent methods are algorithms based on backward induction, partial differential equation 
j—i methods and integral equation methods, cf. e.g. |Det06l Chapter 8]. But these techniques are lim- 

'-H ited to low dimensional problems. Most techniques used today for more complex options are based 

^ on Monte Carlo simulation techniques that were developed in the last years, see JGla04, Chapter 

8] for an overview. We only want to mention stochastic mesh- and regression based-methods, that 
are often combined with using a duality method. 
A not that popular class of algorithms uses the linear programming approach. The basic idea is to 
approximate the underlying process by a Markov chain with a finite state space and to rewrite the 
resulting optimal stopping problem into a linear program and to solve this problem using standard 
techniques, cf. e.g. |CS02] and the references therein for infinite time horizon problems. By the 
curse of dimensionality this approach is limited to low dimensional problems. 
Our approach is different in nature to all approaches described above. The result of this algorithm 
is an analytic approximation to the value function in the continuation set without using discretiza- 
tion. This basic idea is described in the following section and it is shown how semi-infinite linear 
problems come into play. In Section [3] we discuss how the well-known cutting plane algorithm 
can be used to solve such problems. In Section [4] we motivate the further steps by discussing 
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optimal stopping problems with infinite time horizon for one-dimensional diffusion processes. In 
Section [5] we give a theoretical explanation for the accuracy of the algorithm locally around the 
optimization point. The ideas described so fare are then applied to the more interesting case of 
American options on one and more assets with finite time horizon in Section |6] In this section the 
calculation of the Greeks and implied- volatility problems are also discussed. In Section [7] we show 
how the algorithm can be applied for multidimensional diffusions and Levy processes with infinite 
time horizon. Summarizing the results we can say that our algorithm provides good upper bounds 
for the value. In Section |8] we shortly discuss how it can be used to also obtain lower bounds. 
Finally we give a short conclusion in Section [9j 

2 The approach 

We consider a Markovian problem of optimal stopping as follows: 

Let {Xt)t>o be a continuous time strong Markov process with state-space E, g : E ^ [0,oo) 
measurable, T £ (0,oo] and r > 0. We would like to maximize the expected value 

E.,{e-'- g{Xr)) 

over all stopping times t <T with respect to the underlying filtration for all starting points x e E. 
If r = cxD we say that we have an infinite time horizon, if T < cx) we speak of a finite time horizon. 
For convenience we assume T = oo in this section, then the value function does not depend on t; 
this is no real restriction, see |PS06I Chapter I]; in this reference all the following basic facts can 
also be found. 

We define the value function w : i? — > M by 

v{x) = swpE^{e^'"'g{Xr)). 

T 

If we know the value function v, then the optimal stopping problem is solved, but in most situations 
of interest it is not possible to give an explicit expression for v. From the theory for optimal 
stopping Markovian problems it is well-known that under minimal conditions the function v can 
be characterized as the smallest r-excessive function w.r.t. X that majorizes g, i.e. for a fixed 
X{) £ E it holds that 

v(xq) — ini{h{xQ) : h is r-excessive, h > g}. 

Here r-excessive functions are the class of functions, that correspond to the standard supermartin- 
gales for Markov processes, see e.g. |PS06j and the references therein. This formulation corre- 
sponds to the characterization of the value process as the smallest supermartingale dominating the 
gain process in the general setting. Unfortunately the space of r-excessive functions for a process 
X is very wide in general, so that this characterization can be used for an explicit determination 
of the value only in some very special settings. 

Nonetheless using standard terms of optimization theory v{xo) can be seen as the solution to the 
following problem: 

min! h{xo) 

s.t. h{x) > g{x) for all x E E, 
h is r-excessive. 

To see the direct connection to linear programming, let us rewrite the problem as follows: By 
Martin boundary theory - cf. |KW65] - under weak conditions on the process {Xt)t>o each 
r-excessive function h can be represented as 

hi-) = f hi-Mdb), (1) 
Jb 
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where _B is a compact space (called minimal Martin boundary), kb,b g B, are the minimal r- 
excessive functions and tt is a measure on B. Therefore we see that the optimization problem 
described above can be seen as an linear infinite optimization problem: 

h{xo)7r{db) (LIP) 



B 

S.t. / kt,{x)'K{db) > g{x) for all x £ E, 
J B 

TT is a measure on A. 



Here infinite means that we have infinitely many restrictions (since we always assume E to be 
infinite) and we optimize over an infinite-dimensional space of measures. One standard way to 
treat this problem is to discretize the state space E. This leads to an ordinary linear programming 
problem, but this problem is solvable only for low-dimensional spaces E by the curse of dimen- 
sionality. 

The basic idea of our approach is to approximate the value of this linear infinite programming 
problem by reducing the problem to a semi-infinite linear programming problem by choosing a 
finite dimensional subspace of the measure space: 

1. Fix n e N and choose (in a suitable way) a finite subset H := {/ii, /i„} of r-excessive 
functions (equivalently choose n measures tti , . . . , 7r„ ) . 

2. Solve the linear semi-infinite programming problem 

n 

min! ^\,h,{xo) (LSIP) 

1=1 

n 

S.t \ihi{x) > g{x) for all x G E 

1=1 



Remark 2.1. One way for choosing a suitable set H is the following, that will be used in the 
examples below: 

Take a subset H' of the set of all r-excessive functions, that can be parametrized as H' = {ha : a € 
A}. Then choose random parameters a^^-*, a*-"-* G A (e.g. randomly with respect to a suitable 
probability distribution on A) and write hi := ha(i), « = 1, ...,n. 

One immediately obtains the following fact. 

Proposition 2.2. // (A|,...,A*) is a solution of the linear semi-infinite programming problem 



(LSIP), then h* := Y^^=i ^i^ii^) ^•^ upper bound for v{x) for all x E E. 



Proof. The function X]"=i ^l^i r-excessive function majorizing g. By the characterization of 
the value function as smallest r-excessive function majorizing g we obtain the result. □ 

Note that although we only considered one special point xq for the optimization the function 
h* is is an upper bound for the value function v on the whole domain E. As we will see later in 
many situations this is even a good upper bound on a huge neighborhod of xq. But before we 
can apply this algorithm, the first question that arises is how linear semi-infinite programming 



problems of the type ( LSIP I can be solved 



3 Cutting plane method for solving linear semi-infinite pro- 
gramming problems 



The theory for solving linear semi-infinite programming problems of the type ( LSIP I is well de- 



veloped. A good overview is given in [HK93j where theory and algorithms are discussed. One 
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of the key solution techniques is the so-called "cutting plane algorithm". It is based on solving a 
sequence of ordinary linear programming problems, where in each step one further constraint is 
added based on the results obtained so far. To be more precise in step k one considers {xi, Xk} 
and solves 



min! ^X^hiixo) {LPk) 

i=\ 
n 

s.t. \ihi{xj) > g{xj) for each j e {1, fc} 



If \ aIi'''') is an optimal solution to (LPk ) one chooses x^^i as a minimizer of the function 



^Y.^f\,{x)-g{x) 



X 

i=l 



and uses the set {xi, ...,Xk+i} in the next step. This sequence converges to the optimal solution 
under mild restrictions. This algorithm and a relaxed version are e.g. discussed in |WFL98] 

We see that a semi-infinite linear programming problem is reduced to a sequence of standard 
linear programs. The background for the algorithm to work is the reduction theorem for semi- 
infinite linear programs, that states that the infinite restriction set E can be reduced to a set of not 
more than n points, cf. |Kos9H Chapter 10.2]. The points xi,X2,.-. are approximations to these 
points. Therefore the number of iteration steps in the cutting plane algorithm depends on the 
number n of chosen excessive functions and not on the dimension of the underlying space. This is 
the key point for the applicability of our algorithm in higher dimensions. The standard approach 
to solve optimal stopping problems using linear programming is based on the discretization of the 
state space E, cf. e.g. |CS02j . Due to the curse of dimensionality this approach is limited to low 
dimensions. 



For our approach to work we have to find suitable choices for A and H' (see Remark 2.1 1. As 
a motivation for the following sections we first consider one-dimensional diffusion processes with 
infinite time horizon: 



4 One-dimensional diffusions with infinite time horizon 

One-dimensional diffusions have a wide range of applications e.g. in mathematical finance, mathe- 
matical biology, stochastic control and economics. We follow the definition given in |RY99[ Chapter 
VII. 3] that is based on the work of Feller and Ito and McKean (cf. |IM74] ) . i.e. we assume that 
the process is a strong Markov process with continuous sample paths on an interval E. To prevent 
that the interval E can be decomposed into disjoint subintervals from which {Xt)t>a cannot exit, 
we always assume that all diffusions are regular, that is 

Px{Xt — y for some t > 0) > for all x G int{I), y G I. 

To find a suitable set H of r-excessive functions we consider the functions 

^ (x) = /^-(^"''""^{-"<->})' ^ - " 
\ [K(e-'^-=l{,^<^})]-i, x>a 



and 

V'2(a;) 



[^a(e-'^"^l{r.<oo}r\ x<a 
i;a;(e"''^°l{T-„<oo}), x>a, 
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for a fixed point a G int(E). These functions are called the minimal r-harmonic functions. Ob- 
viously ipi is increasing and 'ijj2 is decreasing. Furthermore they are positive, continuous and can 
be used to characterize the boundary behavior of {Xt)t>o. For results in this direction we refer to 
[IM74[ Section 4.6]. All other positive r-harmonic functions are linear combinations of ^pl and ip2- 
In this section we want to study optimal stopping problems for one-dimensional diffusions with 
infinite time horizon. These problems can be solved analytically using different techniques; we 
only refer to |Muc79] . |Sal85] . |BLOO| . |DK03] and |CI10| . Nonetheless in many situations it can 
be helpful to use numerical methods. The following theorem guarantees that H := {1^1,1^2} is a 
reasonable choice for our algorithm to work well. 

Theorem 4.1. Fix xq G E. 

(a) v{xa) is the value of the problem 

min! Ai7/'i(xo) + A2V'2(a;o) 

s.t Xiipi{x) + X2ip2{x) > g{x) for all x d E 



(h) If (Ai, A2) is a solution to the problem above, then v{x) — \i'il)i{x) + \2'4'2{x) for all x in the 
connection component of Xq in the continuation set. 

(c) If an optimal stopping time exists, then 

2 

r, = ini{t > : g{Xt) = ^ X^MXt)} 

i=l 

is optimal under P^. 

Proof, (a) Two proofs based on different methods can be found in |HS10[ Theorem 4.2] and in 
[CiTOl Corollary 2.2]. 



(b) This is a general fact, see Proposition 5.1 in the following section 



(c) If an optimal stopping time exists, then by the general theory the smallest is given by r* — 
inf{< > : v{Xt) = g{Xt)}. Since {Xt)t>o has continuous sample paths the assertion holds by 
the previous point. 

□ 

We consider H := {ipi,ip2}- The theorem states that the value of the linear semi-infinite 



programming problem (LSIPl is equal to v{xo) (and not only an upper bound). Our algorithm 
provides an accurate way for solving these problems using our approach. It is very easy to imple- 
ment in every common language. This and all following examples were implemented in Matlab 



on a standard PC with 1.3 GHz. We used the cutting plane method to solve (LSIPl. This works 
fine and gives the results after some steps of iteration. 

As an example we consider the gain function g(x) = x^ for a standard Brownian motion X. Using 
our algorithm after 5 steps of iteration the linear semi-infinite programming problem reaches 
the solution i;(0) = 5.322 and one furthermore obtains v{x) = 2.661V'i(a:^) + 2.661V'2(a;) for 
X £ [—4.618,4.618], where ipi{x) = e^ '*^''^, "02(2^) — e~°'^'*^^. Moreover the optimal stopping 
time is mi{t > : Xt ^ [—4.618,4.618]}. A graphical illustration can be found in the following 
figure. 

This example is of course not that impressing since optimal stopping problems of this type 
can - in many cases - even be solved analytically by standard techniques such as a free boundary 
approach. But nonetheless this example is instructive for dealing with other problems. We can 
summarize the results as follows: 

• In the definition of the set H one can restrict oneself to r-harmonic functions (instead of 
general r-excessive ones). 
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Figure 1: Example for one-dimensional diffusions and infinite time horizon 



• Optimizing for one point xq in the continuation set yields the value function for the whole 
connection component of the continuation set containing xq- 

• One also obtains the optimal stopping time. 

5 Approximation in the connection component of the con- 
tinuation set 

Next we give the theoretical background for the observation that using an r-harmonic function h 
as an approximation to the value function in a fixed point xq yields a good approximation of the 
value function on the connection component of the continuation set containing xq for a wide class 



of Markov processes (see Theorem 4.1 ^b) above and the numerical results in the next sections): 
Let C denote the connection component of the continuation set that contains xq. We assume that 
the underlying Markov process fulfills a strong maximum principle on C, i.e. we assume that each 
function that is r-harmonic and attains its non-negative maximum in C is constant. This principle 
is well known for certain processes as diffusions under mild conditions, cf. e.g. |Pin95l p. 84], and 
for further classes of processes. Under these assumption we have: 

Proposition 5.1. Let h be an r-harmonic majorant of v with h(xo) = v(xq). 
Then h\c = v\c- 

Proof. Since v is r-harmonic in C so is /i = /i — w. Furthermore h > and h(xo) = 0. Therefore 
by the maximum principle h is constant in C, i.e. h\c = v\c- D 



Now we apply the ideas obtained so fare to the more interesting situation of finite time horizon 
and multidimensional diffusions: 



X 


v{x) 


RLP 


|«(xo)-RLP| 


time 


80 


21.606 


21.615 


0.009 




90 


14.919 


14.923 


0.004 




100 


9.946 


9.951 


0.005 


9.4s 


110 


6.435 


6.439 


0.004 




120 


4.061 


4.064 


0.003 





Table 1: Approximation of the "true" value v{x) (taken form |AC97j ) for the American put 
problem on one asset with time horizon T = 0.5. RLP denotes the values for our algorithms. We 
appHed the algorithm for the optimization point xq = 100 and obtained the other values from this 
optimization as described above. The parameters are d= 1,K = 100, r = 0.06, T = 0.5, a = 0.4. 



6 Time-dependent gain 

In many applications the gain also depends on time, e.g. in mathematical finance one often consid- 
ers problems with a finite time horizon. In this case the value function is an space-time r-excessive 
functions. Before we come to applications we first discuss how the transition densities come into 
play. To this end we use the integral representation as given in ([T]) for this case: 
For a standard Brownian motion it is well known that each excessive function can be written as 
an integral taken over the densities of the Gaussian semigroup (cf . |Sie68j ) . More recently this 
result was extended to a much more general setting in |Jan06] . Different results are given there. 
We only state the following special fact, that is useful for us: 

Under some mild technical assumptions if the underlying transition semigroup is a convolution 
semigroup on a locally compact Abelian group E, then each each space-time excessive function h 
has the representation 



h{t,x) = j l(^s^^){t)pt-s{x,y)'iT{dy,ds), x € E,t>0, 

where tt is a measure on E x [0,oo) and pt{x,y) is a suitably chosen density of the semigroup. 
Standard examples for densities, that are often used are the standard d-dimensional Brownian 
motion, where 

Pt{x,y)^ 
and for Cauchy processes, where 




2t 



With this theoretic result in mind we can treat the well known examples from mathematical 
finance: 



6.1 American put in the Black-Scholes model 

As an example we consider a market Black-Scholes-market where the asset price process X is a 
geometric Brownian motion under the risk neutral measure, that solves 

dXt = rXtdt + aXtdWt 

for a Brownian motion W. Although our approach is also applicable for other gain functions, 
in this subsection we concentrate on the fair price for an American put on X with strike K and 
maturity T as given by 

vit,x)= sup i?(,,,)(e-'^-(if-Xt+,)+), 

T<T-t 
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Figure 2: Graph of the approximated value function h* 



since this example is well studied from different points of view. No closed form solutions are known 
for this problem, but many numerical methods are developed, cf. e.g. |Det06j for an overview. 
The transition density p is given by 

where v = iijcP' ~ 1/2, cf. |BS02] . p. 132. The first idea to apply our algorithms is now to take 
these densities. But one sees that p has a singularity for t = Q. Therefore these densities are 
no good choice, since linear combinations cannot dominate the gain function. Therefore we take 
integrated versions of the density. The easiest such functions are the r-harmonic function given 

by 

/i,(t,a;) = i?(i,,)(e-'-(^-*)l{x^>,}) 

= g-r(T-t)^ / \og{xla) + {r-<j^/2){T-t) \ 
V a^/T^t j 

for a e (0, oo) =: A. Using these functions we can apply our algorithm and compare the results to 
prices taken from |AC97] . The results can be found in the following table. We obtained the data 
in the following way: 

We use our approach with the starting value xq — 100, to = and n ~ 100 and choose the 
parameters ai, ....,a„ according to a uniform distribution on [0, 100]. Applying the optimization 
takes around 10 second. Then we obtained the value at Xq — 100 and tg = and obtained 
parameters AJ^,...,A* such that the function h* :— X]r=i '^i upper bound of the value 
function v. Then we used h* to get the upper bounds for other starting values by just evaluating 
this function at the desired point {t,x). Although h* is optimized for the point {to,xo) = (0, 100) 
the upper bounds for the other points are very good too for other starting prices as shown in 
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X 


interval m [KogOi 




Comp. time 


(80,80) 


[38.01 , 38.35] 


38.30 




(80, 100) 


[32.23 , 32.60] 


32.28 




(80, 120) 


[28.54 , 29.01] 


28.58 




(100,80) 


[33.34 , 33.59] 


33.53 




(100,100) 


[25.81 , 26.02] 


25.86 


42s 


(100,120) 


[20.75 , 21.05] 


20.73 




(120,80) 


[31.21 , 31.31] 


31.30 




(120,100) 


[22.77 , 22.83] 


22.80 




(120,120) 


[16.98 , 16.98] 


16.99 





Table 2: This table states the results for the min-put problem in a Black- Scholes market with two 
assets and parameters K = 100,?' = 0.06, T = 0.5, cti = 0.4, cr2 = 0.8, n — 150. x denotes the 
starting value. We applied our algorithm with starting vale (100,100) and obtained the further 
approximations reported in column "RLP" by evaluating the approximation h* . For comparison 
we also state the values from |Rog02| . 



the table. Hence we only need one approximation for all time horizons and starting values in the 
connected component of the continuation set, see Figure [6?T| Hence we have found an analytically 
given function that is a good approximation to the value function on a huge subset of the time- 
space domain. With the results of Section [5] in mind this is not surprising. 

6.2 American min put on d assets 

As discussed in the introduction it is much more challenging to consider multiple underlyings, i.e. 
the case that X = {X^^\ X^'^^) is a diffusion in a subset of M''. As an example we consider the 
multi-dimensional Black-Scholes market, i.e. X'^^ are geometric Brownian motions with 

fixed correlations of the underlying Brownian motions. One benchmark example in the literature 
is a put option on the minimum of the assets in a Black-Scholes market, i.e. 



g{xi,...,Xd)^\K- min x, 

ie{l,...,d} 



We compare our results to the numerical results given in Rog02 Section 4.2]. With the same 
motivation as for one underlying we could choose the set H' of r-harmonic functions to consist of 
the functions 



where a — (oi, ...,«£() G (OjOo)'*. For highly correlated component processes and high dimensions 
the evaluation of these expectations takes much computational time. In these cases it is more 
reasonable to take the prices of European exchange options between the assets, since these integrals 
can be computed explicitly. Let us remark that all reasonable choices we tried lead to good 
results. We again used the point xq = (100, 100) as the optimization point. Exact results and 
computational time can be found in the following tables. Summarizing we can say 



Optimizing for one special starting point gives very accurate approximations of the value 
function in the continuation set. 



• The same is true for varying time horizons. 



The algorithm also works for large dimensions (e.g. > 10), where normally only Monte-Carlo 
methods are applicable. 
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d 


interval in Rog02 


RLP 


Comp. time 


2 


[24.87 , 25.16] 


24.93 


41 s 


3 


[31.21 , 31.76] 


31.41 


72 s 


4 


[35,72 , 36.28] 


36.01 


115 s 


5 


[39.01 , 39.47] 


39.21 


103 s 


10 


[47.99 , 48.33] 


48.01 


324 s 


15 


[52.23 , 52.14] 


52.10 


612 s 



Table 3: In this table we state the results in the same setting as in Table [6?2| for dimensions 2 to 
15. The parameters now are cr,j — 0.6, T — 0.5, K — 100, r — 0.06 



6.3 Exercise boundary 

Using our approximation of the value function we can also approximate the optimal exercise 
boundary: For the true value function v and each i e [0, T] the exercise boundary is characterized 
as the largest of v(t,-) — g(-). Therefore for an approximation we use the minimizer b(t) of 
{h* — g){t, ■) for each t G [0,T]. A priori it is not clear that this approximation will be accurate 
even if h* is a good approximation of the value function, since 5(-) is very sensitive to the shape of 
h* in both variables x and t. Nonetheless this approximation is very good as shown in the figure 
below. There we compared the boundary to the approximations by standard methods given in 
[LSlOj . We would like to underline that the approximations is even good for small time horizons. 

6.4 Calculation of the Greeks 

For risk management and hedging the Greeks (sensitives) of the option play a major role. The 
Delta of the option - i.e. the derivative with respect to the asset price - is of special interest. 
Using Monte-Carlo techniques it is not straightforward to calculate it, see |Gla04i Chapter 7]. 
Using our method we obtain an approximation h*(t,x) to the value function on the continuation 
set as a function of space and time. Therefore we can calculate the Delta and the Theta of the 
option by simply taking derivatives of h* with respect to x resp. t. Comparing with the results in 
the literature yields that these estimates are very accurate. 

6.5 Calculation of implied volatilities 

Another important topic in the valuation of American options is finding the implied volatilities for 
a given market price vq. From a first view our approach does not seem to be reasonable for this 
question, because the value function for one special volatility first of all does not give information 
about the values for other volatilities. In the following we discuss how this important topic can 
nonetheless be dealt with: 

For a fixed starting volatility CTi we can approximate the price of the asset using our algorithm by 
finding an approximation /ii(0, x\ cti) and can compare this result with the market price Uq. Note 
that this expression gives an explicit function of cti, but for a ^ ai \i is not clear if hi(Q^x\a) 
is an accurate approximation of the price in the model with volatility a. But nonetheless we can 
solve the equation 

hi{Q,x\a) = vq 

for (7. Denote the solution by (T2- Using our approach again we find a new approximation 
/i2(a;, 0; (T2) that can be used to determine 173 and so on. In a general setting there is no hope 
to prove convergence of this sequence to the implied volatility, but nonetheless in our examples 
one obtains a very accurate approximation after three or four steps of iteration even for starting 
volatilities that are fare away from the correct value. This leads to a very easy to implement 
method. Let us emphasize that there is no theoretical justification for the approach to work, but 
nonetheless it seems to work very well. 
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Figure 3: The black graph is the approximated stopping boundary b(-) for parameters T — 1, K — 
100, r — 0.1 obtained by our algorithm. The red graphs are approximations taken from [LSlOj . 
where the upper one is obtained by the PSOR-method and the lower one is the analytical approx- 
imation given by Zho. For comparison the green line is the optimal stopping boundary for the 
perpetual American put. 



7 Infinite time horizon 

After discussing finite time horizon problems in the previous section now we want to discuss the 
case of an infinite time horizon. For practical questions in financial markets this case is not so 
important; perpetual options are only used as a bound for finite time problems. Nonetheless for 
other applications such as sequential statistics and portfolio optimization numerical solutions are 
of importance. Most other numerical methods cannot be applied to these problems, since a dis- 
cretization of an infinite time horizon would be necessary. Exception are the Forward Improvement 
Iteration algorithm discussed in |Irl06] and the results of j CS02] . 

7.1 Multidimensional diffusions 

In the following we are interested in the case that X = {Xt)t>o is a diffusion process with state 
space E C W^,d > 1. For applying our approach we first have to take a suitable subset of 
the class of r-excessive functions w.r.t. X with a suitable parametrization. As in the case of 
finite time horizon we propose to choose a class of r-harmonic functions on E. In this setting 
we propose to take the minimal r-harmonic functions; that are the extreme points of the set of 
all r-harmonic functions on E and can be characterized using the Martin boundary. There is a 
one-to-one correspondence between minimal r-harmonic functions. 

Next we give an examples of interest in mathematical finance where explicit results can be obtained. 
Proposition 7.1. Let X be a d- dimensional Brownian motion on M'' with covariance matrix (aij) 
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and drift ^ = (/ii, ..,^d), i-^- the generator of X is given by 



d 



Write 

^liUi - r = 



d d 

a e E"* : - ^ aijaiaj + ^ , 



{a:: i— ^ exp(a • x) : a € A} 
js the set of all minimal positive r-harmonic functions, where • denotes the usual scalar product. 

This result is well-known, a discussion for this situation can be found in |CI10| . As a standard 
example for a multidimensional problem we consider a perpetual American put option on on index, 
i.e. on a linear combination of assets. This means we consider d (correlated) Brownian motions 

X^^\ X^'^^ with drifts fii,...,iid and covariance structure {(Jij)^ j^i- We interpret e^' ',...,e"^* ' 
as d assets in a Black-Scholes market. Our gain function is given by 

d 

g{xi,...,Xd) = (-FC-^a.e^"')"^ for all {xi,...,Xd) e E"*. 

4 = 1 

Here ai, are positive weight parameters. These options were considered from different points 
of view, see e.g. |Pau01| . 

To use our approach Proposition |7. 1| suggests to take 

A { a e E"* : - ^ aijaiUj + ^ /i^a^ - r = 

4,J = 1 4 = 1 

and 

H' := {x h-> exp(a • x) : a € A}. 

The set A is an ellipsoid and we can choose the random parameters from a uniform 

distribution on A. 

One cannot expect to obtain explicit results for this problem, so that we have to compare our 
results to other approximative results. For this reason we use the forward improvement iteration 
algorithm discussed in ^Irl06j. This algorithm can be applied easily in dimension d = 2, so that 
we use it for our comparison there. We used the forward improvement iteration algorithm with 
a discretization of [0,2] x [0,2] in 100 x 100 points. The results are given in the following table. 
Here the approximation of v{x) using our approach is denoted by RLP, the results by the forward 
improvement iteration algorithm is denoted by FII. 

Although the forward iteration improvement algorithm is limited to low dimensions, our algo- 
rithm is not. It is no problem to to apply it to high dimensional problems. 



7.2 Levy processes with infinite time horizon 

Levy processes are an important class of jump processes that can be used in many fields of 
application such as insurance and finance. Optimal stopping problems with infinite time horizon 
involving Levy processes were studied from different points of view in the last years. For these 
problems overshoot plays a fundamental role. This leads to certain problems for an explicit 
solution. For certain gain functions - such as power functions and functions of put-/call-type - 
semi-explicit solutions were obtained in the terms of the running maximum resp. minimum of the 
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9{xo) 


FII 


RLP 


stopping point 


|RLP-FII| 


(0.7,0.2) 


6.764 


6.764 


6.764 


yes 





(0.7,0.7) 


5.972 


5.976 


5.977 


yes 


0.005 


(1,1) 


4.563 


4.894 


4.944 


no 


0.05 


(1.4,0.6) 


4.122 


4.790 


4.778 


no 


0.012 



Table 4: Results for infinite time horizon and a put on an index. We used the parameters 
Ml = = OjCTij = ^ijiT = 0.1, X = 10, ai = a2 = l,n = 30. With FII we denote the value 
obtained by the forward improvement iteration algorithm and by RLP the values obtained by our 
approach. The optimization toke around 12 seconds. 



process, cf. |Mor02j . (NSOT] , |MS07j and p09] . 

To use our approach we again use the following potential-theoretic representation of r-excessive 
functions: 

As usual we define the resolvent kernel [/'' by 



U'-{x,A) = E, 



J e ^*l[XteA}dt^ , A measurable, x e 



and we assume that U^{x, •) is absolutely continuous with respect to the Lebesgue measure for all 
x G K''. See |Ber96l Chapter 1.3] for a characterization; all the next facts can also be found in this 
reference. 

In the above situation there exists a unique measurable function /i : M"^ — > [0, oo] such that h{- — x) 
is a Lebesgue density of U^{x, ■) for each a; G K'' and y i— >■ h{—y) is r-excessive. For certain 
processes h can be calculated explicitly. Now we can formulate the important representation 
result in the spirit of equation ([ij : 

Proposition 7.2. Any integrahle r-excessive function w can be represented as 

w{x) = / h{a - x)iT{da), x^W'-, 

JtS.'' 

where ir is a unique finite measure on M'^ . 

Remark 7.3. For r = 0, i.e. for problems without discounting, an analogous result holds, if the 
Levy process is assumed to be transient, see JMSO^ for the case d = 1. 

Now we can use the algorithm described above by taking 

A = M'^, H' = {y^h{a-y):aeA}. 

As an example we consider the Novikov-Shiryaev problem, i.e. we use g{x) = (x"'")'*', 7 > 1 as the 
gain function. This problem was completely solved semi-explicitly in |NS04] and |NS07] . 
To check our numerical approach for Levy processes we would like to compare our numerical 
results to explicit results. To this end we assume X to be a compound Poisson process with drift 
and positive exponential jumps, i.e. X has the form 

Xt = ct + J2Yz, i>0, 

i=l 

where c < 0, {Nt)t>o is a Poisson process with intensity A and (li)igN is a sequence of independent 
i?a;p(a)-distributed random variables. In this setting an explicit solution was obtained in jMS07] . 
In this case the Green function h is given by 



h{x) 



A2eP^, x<Q 
-Ai, x>Q, 
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Figure 4: Gain function (black) and approximated value function (red) in the Levy example for 
parameters r — 2,a ^ 1,P = 0.5, A/c = —0.5, 7 = 2.5. 



where p a + A/c > 0, Ai := a/{X + ca) and A2 := \/{aX + a'^a), see |MS07I Section 5]. 



We applied our algorithm with n ~ 150, and choosing the parameters a' 



(1) 



Jn) 



according to a 



uniform distribution on the interval [0, 20] turned out to be a reasonable choice. The computational 



results are given in Figure 7.2 Using these r-excessive functions we obtain a good approximation 
not only on the continuation set, but also on the optimal stopping set. The valuation of American 
options with finite time horizon in Levy markets will be discussed in detail in a forthcoming paper. 



8 Lower bounds of the value function 



As explained above our method immediately leads to good upper bounds of the value function. 
This is indeed the important contribution of our approach since most Monte-Carlo methods leads 
to good lower bounds, but nonetheless to deal with new problems one also would like to obtain 
lower bounds for the value. In this section we discuss how this can be realized using our approach. 
For easy examples like the American put in the Black-Scholes market the method also gives an 
approximation to the stopping boundary. Using the stopping time associated with this boundary 
leads to very good lower bounds. 

For more complex examples the idea is to use the approximation of the value function h* = 
Sr=i K^i- The first idea is to choose the optimal stopping time 

T* =mf{t>0: g{Xt)>vit,Xt)} 

and to substitute v hy h* . But since h* > g and h* is just an approximation to v this stopping 
time does not seem to be appropriate. Instead we choose e > and take 

r' = mf{t > : g{Xt) + e > v{t, Xt)}, 
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then it is well known that r' is e-optimal in the sense that 

E^t,x){e-'"^'g{Xr')) > v{t,x) - e 
under minimal conditions. Now we take 

TO := inf{t > : g{Xt) + e > h*{t,Xt)} 

and use 

is a lower bound of v{t,x). In most problems we cannot find analytical expressions for the expec- 
tation on the right hand side, but it can be approximated using Monte Carlo techniques. In our 

examples the lower approximations were quite near the upper ones. 

Another approach is to use the r-harmonic function h* for variants of other methods like the 
Longstaff-Schwartz algorithm. This will be discussed by the author in a forthcoming paper. 

9 Conclusion 

As a conclusion let us summarize the important properties of the approach described above: 

• The algorithm is based on reducing the ILP-problem connected to optimal stopping to a 
SILP-problem by choosing finitely many r-excessive functions. 

• No discretization of space and time and no Monte-Carlo-elements are necessary. 

• The algorithm is very easy to implement in every common language (1 page of programming 
code!). 

• Optimizing for one special starting point gives very accurate upper bound of the value 
function in the continuation set. 

• The same is true for varying time horizons. 

• The algorithm also works for large dimensions (e.g. > 10), where apart from it only Monte- 
Carlo methods arc applicable. 

• The Greeks can be found immediately. 

• Implicit volatilities can be calculated. 

• The algorithm can be used for infinite time horizons, too. 
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